## 2014 Data

setwd("C:/Users/Nathen Huang/Dropbox/QMSS Thesis/Martha's Vineyard/Data/Ready for Analysis")
mv14 <- read.csv("MV 2014 Data for Analysis.csv")
mv14prac <- read.csv("MV 2014 Data for Analysis.csv")
table(mv14prac$oppolice)
mv14$X <- NULL
mv14$X.1 <- NULL
mv14$year <- rep(2014)

which(colnames(mv14)=="lgencor1")
which(colnames(mv14)=="lpres9")

mv14sub <- mv14[,7:15]
mv14sub <- mv14sub[,-c(5,6)]
colnames(mv14sub)
cronbach(mv14sub)

mv14$legal <- rowMeans(mv14[,7:15], na.rm = TRUE)
summary(mv14$legal)


outcomes14 <- mv14[,7:15]

mv14$lgencor1 <- as.numeric(mv14$lgencor1)
mv14$lsupcor2 <- as.numeric(mv14$lsupcor2)
mv14$lpolice3 <- as.numeric(mv14$lpolice3)
mv14$llocpol4 <- as.numeric(mv14$llocpol4)
mv14$lcong7 <- as.numeric(mv14$lcong7)
mv14$llocgov8 <- as.numeric(mv14$llocgov8)
mv14$lpres9 <- as.numeric(mv14$lpres9)

outcomes14_fa <- factanal(na.omit(outcomes14), factors = 5, rotation = "varimax")
outcomes14_fa

Legal System 2004
-Local Police
-Court System 
-Local Police/President 
-Congress/Local Government 


Legal System 2014
-police/local police
-Supreme Court/Court System 
-Congress
-The President

Legal System Merge
-Supreme Court/Court System/Congress 
-Local Police/National Police
-President 
-Local Government 


mv14$doccup <- as.character(mv14$doccup)
str(mv14$doccup)

tail <- tail(merge)

table(merged2$year)

mv14$demracecheck <- as.numeric(mv14$demracecheck) 
mv14$race.check <- ifelse(mv14$demracecheck == 1, 1, ifelse(mv14$demracecheck == 2 | mv14$demracecheck == 3, 2, 3))
mv14$race.check <- as.factor(mv14$race.check)
table(mv14$race.check)

## Selecting white respondents (using very specific criteria)
mv14$white <- ifelse(mv14$demracecheck == 1, 1, 0)
mv14$white.check <- ifelse(mv14$demracecheck == 1, 1, 0)

## Selecting black respondents (Only if Black respondents self-reported Black identity)
mv14$black <- ifelse(mv14$demracecheck == 2 | mv14$demracecheck == 3, 1, 0)
mv14$black.check <- ifelse(mv14$demracecheck == 2 | mv14$demracecheck == 3, 1, 0)

table(mv14$white.check) ## 331 White respondents 

table(mv14$black.check) ## 232 Black respondents 

mv14$poor <- ifelse(mv14$dincome < 4, 1, 0)
mv14$rich <- ifelse(mv14$dincome > 6, 1, 0)

white.black.only.14 <- filter(mv14, white.check == 1 | black.check == 1)
white.black.only.14$white.black <- ifelse(white.black.only.14$white.check == 1, "white", "black")
white.black.only.14$poor.rich <- ifelse(white.black.only.14$poor == 1, "poor", ifelse(white.black.only.14$rich == 1, "rich", NA))

lpolicec <- summarySE(white.black.only.14, measurevar="lpolice3", groupvars=c("white.black","poor.rich"), na.rm = TRUE)
lpolicec <- lpolicec[-c(3,6),]
lpresc$row.names <- NULL

white.black.only.14$white.black <-as.factor(white.black.only.14$white.black)
white.black.only.14$poor.rich <-as.factor(white.black.only.14$poor.rich)

# Error bars represent standard error of the mean
ggplot(lpolicec, aes(x=poor.rich, y=lpolice3, fill=white.black)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=lpolice3-se, ymax=lpolice3+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9)) +
  scale_fill_manual(values = c("black","bisque","black","bisque"),
                    name="Race of Respondent",
                    breaks=c("black", "white"),
                    labels=c("Black", "White")) +
  xlab("Class Categories") +
  ylab("Favorability of National Police") +
  ggtitle("Favorability of National Police, 2014") +
  scale_y_continuous(breaks=seq(0, 3.5, 0.5)) +
  theme_bw()

llocpolc <- summarySE(white.black.only.14, measurevar="llocpol4", groupvars=c("white.black","poor.rich"), na.rm = TRUE)
llocpolc <- llocpolc[-c(3,6),]
lpresc$row.names <- NULL

white.black.only.14$white.black <-as.factor(white.black.only.14$white.black)
white.black.only.14$poor.rich <-as.factor(white.black.only.14$poor.rich)

# Error bars represent standard error of the mean
ggplot(llocpolc, aes(x=poor.rich, y=llocpol4, fill=white.black)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=llocpol4-se, ymax=llocpol4+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9)) +
  scale_fill_manual(values = c("black","bisque","black","bisque"),
                    name="Race of Respondent",
                    breaks=c("black", "white"),
                    labels=c("Black", "White")) +
  xlab("Class Categories") +
  ylab("Favorability of Local Police") +
  ggtitle("Favorability of Local Police, 2014") +
  scale_y_continuous(breaks=seq(0, 3.5, 0.5)) +
  theme_bw()

mv14$white.black <- ifelse(mv14$white.check == 1, 1, 
                           ifelse(mv14$black.check == 1, 2, NA))

mv14$poor.rich <- ifelse(mv14$dincome < 4, 1, 
                         ifelse(mv14$dincome > 6, 2, NA))

table(mv14$poor.rich)


## DO NOT DELETE THIS CODING SCHEME OMG THIS TOOK FOREVER
mv14$type <- ifelse(mv14$white.check == 1, ifelse(mv14$poor == 1, "poor.white", ifelse(mv14$rich == 1, "rich.white", NA)), ## TRUE STATEMENT
                    ifelse(mv14$black.check == 1, ifelse(mv14$poor == 1, "poor.black", ifelse(mv14$rich == 1, "rich.black", NA)), NA))

table(mv14$type)


## Perceived Opportunity 
mv14$opeq <- 5 - mv14$opeq ## Reverse-code- higher values = higher perception of fairness


opp.disc <- c("ophouse", "oppolice")
table(mv14$oppolice)
mv14$ophouse <- 5 - as.numeric(mv14$ophouse)
summary(mv14$ophouse)
mv14$oppolice <- 6 - as.numeric(mv14$oppolice)
summary(mv14$oppolice)
mv14$op.disc <- rowMeans(mv14[,opp.disc], na.rm = TRUE)
mv04$op.disc <- rowMeans(mv04[,opp.disc], na.rm = TRUE)

mv14$dempolt <- as.factor(mv14$dempolt)
mv14$demskin <- as.numeric(mv14$demskin)
mv14$demgend <- as.factor(mv14$demgend)

mv14$edcourtesy <- 7 - mv14$edcourtesy 
mv14$edinsult <- 7 - mv14$edinsult
mv14$edservc <- 7 - mv14$edservc 
mv14$edbettr <- 7 - mv14$edbettr  
mv14$ednsmart <- 7 - mv14$ednsmart

every.disc <- c("edcourtesy", "edinsult", "edservc", "edbettr", "ednsmart")
mv14$mean.disc <- rowMeans(mv14[,every.disc], na.rm = TRUE)

## Visualizing differences

bwplot(type~police.trust,data=mv14,
       xlab="Mean Perceived Favorability of Police",
       main="Favorable Views of Police Given Race/Class, 2014")

bwplot(type~mean.govtrust,data=mv14,
       xlab="Mean Perceived Favorability of Legal System",
       main="Favorable Views of Legal System Given Race/Class, 2014")

bwplot(type~opeq,data=mv14,
       xlab="Mean Perceived Opportunity in the US",
       main="Perception of Opportunity in the US, 2014")

y <- c(mv04$mean.rrs, mv14$mean.disc)

group <- as.factor(c(rep(1, length(mv04$mean.rrs)), rep(2, length(mv14$mean.disc))))

leveneTest(y, group)

outcomes14 <- na.omit(blacks.14[,dependents])
str(outcomes14)

outcomes14$lpubsch5 <- as.numeric(outcomes14$lpubsch5)
outcomes14$llocsch6 <- as.numeric(outcomes14$llocsch6)

outcomes_fa_14 <- factanal(outcomes14, factors = 5, rotation = "varimax")

outcomes_fa_14

## Dividing population by Black and White respondents 
library(dplyr)


blacks.14 <- filter(mv14, black.check == 1)
whites.14 <- filter(mv14, white.check == 1)

summary(mv14$op.disc)

t.test(whites.14$op.disc, blacks.14$op.disc)
cohen.d(whites.14$op.disc, blacks.14$op.disc, na.rm = TRUE)

t.test(blacks.14$legal, whites.14$legal)
cohen.d(blacks.14$lega, whites.14$legal, na.rm = TRUE)

table(mv14$type)

t.test(whites.14$op.disc, blacks.14$op.disc)
cohen.d(whites.14$op.disc, blacks.14$op.disc, na.rm = TRUE)

poor.14 <- filter(mv14, dincome < 4)
rich.14 <- filter(mv14, dincome > 6)

poor.blacks.14 <- filter(blacks.14, dincome < 4)
rich.blacks.14 <- filter(blacks.14, dincome > 6)

t.test(poor.blacks.14$op.disc, rich.blacks.14$op.disc)
str(mv14$op.disc)
summary(mv14$op.disc)
cohen.d(poor.blacks.14$op.disc, rich.blacks.14$op.disc, na.rm = TRUE)

poor.whites.14 <- filter(whites.14, dincome < 4)
rich.whites.14 <- filter(whites.14, dincome > 6)

library(effsize)


library(ggplot2)
par(mfrow=c(2,3))


# Plots
library(ggplot2)
colnames(mv14)
dfwc <- summarySEwithin(mv04, measurevar="value", withinvars="condition",
                        subject.id="subject", na.rm=FALSE, conf.interval=.95)
install.packages("Rmisc")
library(Rmisc)
tgc <- summarySE(mv14$lpres9)

which(colnames(mv14) == "lpres9")
mv14sub <- mv14[,7:15]

str(mv14$lpres9)
mv14$lpresSE <- 1.96*(sd(mv14$lpres9, na.rm=T)/sqrt(661))

tgc <- summarySE(white.black.only, measurevar="lpres9", groupvars=c("white.black","poor.rich"), na.rm = TRUE)


df2 <- data_summary(white.black.only, varname = "lpres9", groupnames = c("white.black", "poor.rich"))
df2$white.black <- as.factor(df2$white.black)
df2$poor.rich <- as.factor(df2$poor.rich)
  
summary(mv14$lpres9)

ggplot(mv14[, c("type", "lpres9")], 
       aes(x=type, y=lpres9, fill = type)) +  
  geom_bar() +
  scale_fill_manual(values = c("black","bisque","black","bisque")) +
  ggtitle("Favorability of National Police, 2014") +
  theme_bw()


ggplot(na.omit(mv14[, c("type", "lpolice3")]), 
       aes(x=type, y=lpres9, fill = type)) +  
  geom_bar() +
  geom_errorbar(aes(ymin=lpres9-lpresSE, ymax=lpolice3+lpresSE),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9)) +
  scale_fill_manual(values = c("black","bisque","black","bisque")) +
  stat_summary(fun.y="mean", geom="bar") +
  ggtitle("Favorability of National Police, 2014") +
  theme_bw()

ggplot(na.omit(mv14[, c("type", "llocpol4")]), 
       aes(x=type, y=llocpol4, fill = type)) +  
  geom_bar() +
  scale_fill_manual(values = c("black","bisque","black","bisque")) +
  stat_summary(fun.y="mean", geom="bar") +
  ggtitle("Favorability of Local Police, 2014") +
  theme_bw()

ggplot(na.omit(mv14[, c("type", "lpres9")]), 
       aes(x=type, y=lpres9, fill = type)) +  
  geom_bar() +
  geom_errorbar(aes(ymin=df$x-df$se, ymax=df$x+df$se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9)) +
  scale_fill_manual(values = c("black","bisque","black","bisque")) +
  stat_summary(fun.y="mean", geom="bar") +
  ggtitle("Favorability of Local Police, 2014") +
  theme_bw()

lpresc <- summarySE(white.black.only, measurevar="lpres9", groupvars=c("white.black","poor.rich"), na.rm = TRUE)
lpresc <- lpresc[-c(3,6),]
lpresc$row.names <- NULL

white.black.only$white.black <-as.factor(white.black.only$white.black)
white.black.only$poor.rich <-as.factor(white.black.only$poor.rich)

# Error bars represent standard error of the mean
ggplot(lpresc, aes(x=poor.rich, y=lpres9, fill=white.black)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=lpres9-se, ymax=lpres9+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9)) +
  scale_fill_manual(values = c("black","bisque","black","bisque"),
                    name="Race of Respondent",
                    breaks=c("black", "white"),
                    labels=c("Black", "White")) +
  xlab("Class Categories") +
  ylab("Favorability of the President") +
  ggtitle("Race and Wealth Differentiation in Perceptions of President Favorability, 2004")

bp + 
  
ToothGrowth$dose.cat <- factor(ToothGrowth$dose, labels=paste("d", 1:3, sep=""))
df <- with(white.black.only, aggregate(lpres9, list(white.black=white.black, poor.rich=poor.rich), mean, na.rm=TRUE))
df$se <- with(white.black.only , aggregate(lpres9, list(white.black=white.black, poor.rich=poor.rich), 
                                      function(x) sd(na.omit(x))/sqrt(1021)))[,3]

ggplot(df, 
       aes(x=poor.rich, y=x, fill = white.black)) +  
  geom_bar() +
  geom_errorbar(aes(ymin=x-se, ymax=x+se),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9)) +
  scale_fill_manual(values = c("black","bisque","black","bisque")) +
  ggtitle("Favorability of Local Police, 2014") +
  theme_bw()

ggplot(na.omit(mv14[, c("type", "lpres9", "lpresSE")]), 
       aes(x=type, y=lpres9, fill = type)) +  
  geom_bar() +
  geom_errorbar(aes(ymin=lpres9-lpresSE, ymax=lpres9+lpresSE),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9)) +
  scale_fill_manual(values = c("black","bisque","black","bisque")) +
  stat_summary(fun.y="mean", geom="bar") +
  ggtitle("Favorability of Local Police, 2014") +
  theme_bw()

t.test(blacks.04$mean.govtrust, blacks.14$mean.govtrust)
t.test(blacks.04$police.trust, blacks.14$police.trust)
t.test(blacks.04$opeq, blacks.14$opeq)
t.test(blacks.04$dwpimp, blacks.14$dwpimp)

t.test(blacks.04$dwpimp, blacks.14$dwpimp)

table(blacks.04$dwpimp)
table(blacks.14$dwpimp)

par(mfrow=c(1,2))
plot(mean.govtrust ~ white.black + poor.rich, data=mv14)

## 1. lgencor1 

library(effsize)

t.test(poor.14$lgencor1, rich.14$lgencor1)
cohen.d(poor.14$lgencor1, rich.14$lgencor1, na.rm = TRUE)


t.test(whites.14$lgencor1, blacks.14$lgencor1)
sd(whites.14$lgencor1, na.rm = TRUE)
sd(blacks.14$lgencor1, na.rm = TRUE)

cohen.d(whites.14$lgencor1, blacks.14$lgencor1, na.rm = TRUE)
table(mv14$type)

t.test(poor.blacks.14$lgencor1, rich.blacks.14$lgencor1)
sd(poor.blacks.14$lgencor1, na.rm = TRUE)
sd(rich.blacks.14$lgencor1, na.rm = TRUE)
cohen.d(poor.blacks.14$lgencor1, rich.blacks.14$lgencor1, na.rm = TRUE)

t.test(poor.whites.14$lgencor1, rich.whites.14$lgencor1)
cohen.d(poor.whites.14$lgencor1, rich.whites.14$lgencor1, na.rm = TRUE)


results1.14 = lm(lgencor1 ~ white.black*poor.rich, data= white.black.only.14)
aov1.14 <- aov(results1)
anova(results1.14)
aov1.14
TukeyHSD(aov1.14)
anova(na.omit(results1))

mv14$dempolt <- as.factor(mv14$dempolt)
mv14$demskin <- as.numeric(mv14$demskin)
mv14$demgend <- as.factor(mv14$demgend)

ols1 <- lm(lgencor1 ~ op.disc + demskin + race.check*dincome + demage + demgend + dempolt + dedu.10, data = mv14)
summary(ols1)
bptest(ols1)
coeftest(ols1, vocv = hccm(ols1)) ## Robust Standard Errors Correction

## 2. lsupcor2

t.test(poor.14$lsupcor2, rich.14$lsupcor2)
cohen.d(poor.14$lsupcor2, rich.14$lsupcor2)

t.test(whites.14$lsupcor2, blacks.14$lsupcor2)
sd(whites.14$lsupcor2, na.rm = TRUE)
sd(blacks.14$lsupcor2, na.rm = TRUE)

cohen.d(whites.14$lsupcor2, blacks.14$lsupcor2, na.rm = TRUE)
sd(poor.blacks.14$lsupcor2, na.rm = TRUE)
sd(rich.blacks.14$lsupcor2, na.rm = TRUE)
t.test(poor.blacks.14$lsupcor2, rich.blacks.14$lsupcor2)

cohen.d(poor.blacks.14$lsupcor2, rich.blacks.14$lsupcor2, na.rm = TRUE)

t.test(poor.whites.14$lsupcor2, rich.whites.14$lsupcor2)

results2 = lm(lsupcor2 ~ white.black*poor.rich, data= white.black.only.14)
aov2.14 <- aov(results2)
TukeyHSD(aov2.14)
anova(na.omit(results2))

ols2 <- lm(lsupcor2 ~ op.disc + demskin + race.check*dincome + demage + demgend + dempolt + dedu.10, data = mv14)
summary(ols2)
bptest(ols2)
coeftest(ols2, vocv = hccm(ols2)) ## Robust Standard Errors Correction

## 3. lpolice3

t.test(poor.14$lpolice3, rich.14$lpolice3)
cohen.d(poor.14$lpolice3, rich.14$lpolice3)

t.test(whites.14$lpolice3, blacks.14$lpolice3)
sd(whites.14$lpolice3, na.rm = TRUE)
sd(blacks.14$lpolice3, na.rm = TRUE)

cohen.d(whites.14$lpolice3, blacks.14$lpolice3, na.rm = TRUE)

t.test(poor.blacks.14$lpolice3, rich.blacks.14$lpolice3)
sd(poor.blacks.14$lpolice3, na.rm = TRUE)
sd(rich.blacks.14$lpolice3, na.rm = TRUE)
cohen.d(poor.blacks.14$lpolice3, rich.blacks.14$lpolice3, na.rm = TRUE)

t.test(poor.whites.14$lpolice3, rich.whites.14$lpolice3)
sd(poor.whites.14$lpolice3, na.rm = TRUE)
sd(rich.whites.14$lpolice3, na.rm = TRUE)
cohen.d(poor.whites.14$lpolice3, rich.whites.14$lpolice3, na.rm = TRUE)

results3 = lm(lpolice3 ~ white.black * poor.rich, data= white.black.only.14)
aov3.14 <- aov(results3)
TukeyHSD(aov3.14)

ols3 <- lm(lpolice3 ~ op.disc + demskin + race.check*dincome + demage + demgend + dempolt + dedu.10, data = mv14)
summary(ols3)
bptest(ols3)
coeftest(ols3, vocv = hccm(ols3)) ## Robust Standard Errors Correction

# 4. llocpol4 

t.test(poor.14$llocpol4, rich.14$llocpol4)
cohen.d(poor.14$llocpol4, rich.14$llocpol4, na.rm = TRUE)

t.test(whites.14$lcong7, blacks.14$llocpol4)
sd(whites.14$llocpol4, na.rm = TRUE)
sd(blacks.14$llocpol4, na.rm = TRUE)
cohen.d(whites.14$llocpol4, blacks.14$llocpol4, na.rm = TRUE)

t.test(whites.14$op.disc, blacks.14$op.disc)
sd(whites.14$op.disc)
sd(blacks.14$op.disc, na.rm = TRUE)
cohen.d(whites.14$op.disc, blacks.14$op.disc, na.rm = TRUE)

t.test(poor.blacks.14$op.disc, rich.blacks.14$op.disc)
sd(poor.blacks.14$op.disc, na.rm = TRUE)
sd(rich.blacks.14$op.disc, na.rm = TRUE)
cohen.d(poor.blacks.14$op.disc, rich.blacks.14$op.disc, na.rm = TRUE)

t.test(poor.blacks.14$llocpol4, rich.blacks.14$llocpol4)
sd(poor.blacks.14$llocpol4, na.rm = TRUE)
sd(rich.blacks.14$llocpol4, na.rm = TRUE)
cohen.d(poor.blacks.14$llocpol4, rich.blacks.14$llocpol4, na.rm = TRUE)

t.test(poor.blacks.14$llocpol4, rich.blacks.14$llocpol4)
sd(poor.blacks.14$llocpol4, na.rm = TRUE)
sd(rich.blacks.14$llocpol4, na.rm = TRUE)
cohen.d(poor.blacks.14$llocpol4, rich.blacks.14$llocpol4, na.rm = TRUE)

t.test(whites.04$mean.rrs, blacks.04$mean.rrs)
sd(whites.04$mean.rrs, na.rm = TRUE)
sd(blacks.04$mean.rrs, na.rm = TRUE)
cohen.d(whites.04$mean.rrs, blacks.04$mean.rrs, na.rm = TRUE)

t.test(poor.whites.14$llocpol4, rich.whites.14$llocpol4)
sd(poor.whites.14$llocpol4, na.rm = TRUE)
sd(rich.whites.14$llocpol4, na.rm = TRUE)
cohen.d(poor.whites.14$llocpol4, rich.whites.14$llocpol4, na.rm = TRUE)

results4.14 = lm(llocpol4 ~ white.black * poor.rich, data = white.black.only.14)
aov4.14 <- aov(results4.14)
anova(results4.14)
TukeyHSD(aov4.14)
anova(na.omit(results4))




ols4 <- lm(llocpol4 ~ op.disc + demskin + race.check*dincome + demage + demgend + dempolt + dedu.10, data = mv14)
summary(ols4)
bptest(ols4)
coeftest(ols4, vocv = hccm(ols4)) ## Robust Standard Errors Correction


# 5. lcong7

t.test(poor.14$lcong7, rich.14$lcong7)
cohen.d(poor.14$lcong7, rich.14$lcong7)

t.test(whites.14$lcong7, blacks.14$lcong7)
sd(whites.14$lcong7, na.rm = TRUE)
sd(blacks.14$lcong7, na.rm = TRUE)

cohen.d(whites.14$lcong7, blacks.14$lcong7, na.rm = TRUE)

t.test(poor.blacks.14$lcong7, rich.blacks.14$lcong7)
sd(poor.blacks.14$lcong7, na.rm = TRUE)
sd(rich.blacks.14$lcong7, na.rm = TRUE)
cohen.d(poor.blacks.14$lcong7, rich.blacks.14$lcong7, na.rm = TRUE)

t.test(poor.whites.14$lcong7, rich.whites.14$lcong7)

results7 = lm(lcong7 ~ white.black*poor.rich, data= white.black.only.14)
aov7.14 <- aov(results7)
TukeyHSD(aov7.14)
anova(na.omit(results7))

ols7 <- lm(lcong7 ~ op.disc + demskin + race.check*dincome + demage + demgend + dempolt + dedu.10, data = mv14)
summary(ols7)
bptest(ols7)
coeftest(ols7, vocv = hccm(ols7)) ## Robust Standard Errors Correction

# 6. llocgov 8

t.test(poor.14$llocgov8, rich.14$llocgov8)
cohen.d(poor.14$llocgov8, rich.14$llocgov8)

t.test(whites.14$llocgov8, blacks.14$llocgov8)
sd(whites.14$llocgov8, na.rm = TRUE)
sd(blacks.14$llocgov8, na.rm = TRUE)


cohen.d(whites.14$llocgov8, blacks.14$llocgov8, na.rm = TRUE)

t.test(poor.blacks.14$llocgov8, rich.blacks.14$llocgov8)
sd(poor.blacks.14$llocgov8, na.rm = TRUE)
sd(rich.blacks.14$llocgov8, na.rm = TRUE)
cohen.d(poor.blacks.14$llocgov8, rich.blacks.14$llocgov8, na.rm = TRUE)

t.test(poor.whites.14$llocgov8, rich.whites.14$llocgov8)

results8 = lm(llocgov8 ~ white.black * poor.rich, data = white.black.only.14)
aov8.14 <- aov(results8)
TukeyHSD(aov8.14)
anova(na.omit(results8))

ols8 <- lm(llocgov8 ~ op.disc + demskin + race.check*dincome + demage + demgend + dempolt + dedu.10, data = mv14)
summary(ols8)
bptest(ols8)
coeftest(ols8, vocv = hccm(ols8)) ## Robust Standard Errors Correction

# 7. lpres9
library(stargazer)
t.test(poor.14$lpres9, rich.14$lpres9)
cohen.d(poor.14$lpres9, rich.14$lpres9)

stargazer(ols1, title="Results", align=TRUE)
install.packages('dcolunn')

t.test(whites.14$lpres9, blacks.14$lpres9)
sd(whites.14$lpres9, na.rm = TRUE)
sd(blacks.14$lpres9, na.rm = TRUE)
cohen.d(whites.14$lpres9, blacks.14$lpres9, na.rm = TRUE)

t.test(poor.blacks.14$lpres9, rich.blacks.14$lpres9)
sd(poor.blacks.14$lpres9, na.rm = TRUE)
sd(rich.blacks.14$lpres9, na.rm = TRUE)
cohen.d(poor.blacks.14$lpres9, rich.blacks.14$lpres9, na.rm = TRUE)

t.test(poor.whites.14$lpres9, rich.whites.14$lpres9)

results9 = lm(lpres9 ~ white.black*poor.rich, data = white.black.only.14)
aov9.14 <- aov(results9)
TukeyHSD(aov9.14)
anova(na.omit(results9))

results10 = lm(mean.rrs ~ white.black*poor.rich, data = white.black.only)
aov10 <- aov(results10)
TukeyHSD(aov10)

ols9 <- lm(lpres9 ~ op.disc + demskin + race.check*dincome + demage + demgend + dempolt + dedu.10, data = mv14)
summary(ols9)
bptest(ols9)
coeftest(ols9, vocv = hccm(ols9)) ## Robust Standard Errors Correction